Criminalità e immigrazione: fenomeni in relazione o strumento di proaganda?

library(tidyverse)
library(readxl)
library(plotly)
library(countrycode)
library(modelr)
library(broom)
library(ggrepel)

Uno dei temi politici e sociologici più caldo e ricorrente degli ultimi tempi è stato ed è ancora sicuramente il fenomeno dell’immigrazione, non a caso questa tematica è stata anche oggetto di propaganda per alcuni partiti politici. Spesso l’aumento in volume di tale fenomeno è visto con timore e preoccupazione, non solo per il costo di integrazione ma anche perchè in molti ipotizzano che vi sia un qulche tipo di legame fra criminalità e immigrazione. Ma la domanda allora è: esiste davvero un legame diretto fra queste due entità o è solo una percezione fuorviata dai mass media e dalla propaganda politica? Lo scopo dell’analisi è rispondere a questo interrogativo andando ad analizzare congiuntamente database sulla criminalità sull’immigrazione. Per quanto riguarda la criminalità si è deciso di concentrarsi sui seguenti 3 tipi di crimine:

• Omicidio
• Aggressione
• Violenza sessuale

1 Fonti di informazione

• Per la criminalità:
• World Bank
• UNODC
• Per l’immigrazione:
• Nazioni unite

2 Importazione dei dati

Nel seguito si procederà con l’importazione dei file ottenuti dai siti elencati in precedenza.

2.1 Importazione file su omicidi

Procediamo quindi con l’importazione dei dati: Importiamo il dataset sugli omicidi ottenuto da WorldBank

omicidiWB = read_csv("OmicidiWB.csv", na = "", skip = 4)
omicidiWB

# Il dataset non è in formato tidy in quanto abbiamo una colonna per ogni anno di osservzione,
# la cosa da fare è quindi rendere il database nel formto adeguato
unique(omicidiWB$X63)

## [1] NA

omicidiWBT = gather(omicidiWB, Year, Value, 5:62)
glimpse(omicidiWBT)

## Observations: 15,312
## Variables: 7
## $ `Country Name`   <chr> "Aruba", "Afghanistan", "Angola", "Albania", ...
## $ `Country Code`   <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARB", "AR...
## $ `Indicator Name` <chr> "Intentional homicides (per 100,000 people)",...
## $ `Indicator Code` <chr> "VC.IHR.PSRC.P5", "VC.IHR.PSRC.P5", "VC.IHR.P...
## $ X63              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ Year             <chr> "1960", "1960", "1960", "1960", "1960", "1960...
## $ Value            <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...

omicidiWBT = omicidiWBT[,-5]
# i nomi di alcune colonne sono composti da due parole separate da uno spazio, per comodità sostituiamo lo spazio con il carattere "_"
colnames(omicidiWBT) = gsub(" ", "_", colnames(omicidiWBT))
# andiamo a vedere il risultato della correzione
glimpse(omicidiWBT)

## Observations: 15,312
## Variables: 6
## $ Country_Name   <chr> "Aruba", "Afghanistan", "Angola", "Albania", "A...
## $ Country_Code   <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARB", "ARE"...
## $ Indicator_Name <chr> "Intentional homicides (per 100,000 people)", "...
## $ Indicator_Code <chr> "VC.IHR.PSRC.P5", "VC.IHR.PSRC.P5", "VC.IHR.PSR...
## $ Year           <chr> "1960", "1960", "1960", "1960", "1960", "1960",...
## $ Value          <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...

#vediamo che le variabili Year e Value, che dovrebbero esere in formato numerico sono in formato character. procediamo quindi alla conversione
omicidiWBT$Year = as.integer(omicidiWBT$Year)
omicidiWBT$Value = as.numeric(omicidiWBT$Value)

#risultato finale
summary(omicidiWBT)

##  Country_Name       Country_Code       Indicator_Name    
##  Length:15312       Length:15312       Length:15312      
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##                                                          
##  Indicator_Code          Year          Value        
##  Length:15312       Min.   :1960   Min.   :  0.000  
##  Class :character   1st Qu.:1974   1st Qu.:  1.487  
##  Mode  :character   Median :1988   Median :  3.700  
##                     Mean   :1988   Mean   :  8.418  
##                     3rd Qu.:2003   3rd Qu.:  9.500  
##                     Max.   :2017   Max.   :139.132  
##                                    NA's   :12693

head(omicidiWBT)

La struttura finale prevede un dataset di 15576 osservazioni e 6 variabili:

1. Country_Name: il nome del paese
2. Country_Code: il codice identificativo del paese
3. Indicator_name: il nome dell’indicatore cosiderato, ovvero il numero di casi di omicidio su 100000 persone
4. Indicator_Code: codice identificativo dell’indicatore
5. Year: Anno relativo all’osservazione
6. Value: valore dell’indicatore

Andiamo ad analizzare il contenuto delle variabili Indicator_Name e Indicator_code

unique(omicidiWBT$Indicator_Name)

## [1] "Intentional homicides (per 100,000 people)"

unique(omicidiWBT$Indicator_Code)

## [1] "VC.IHR.PSRC.P5"

Come possiamo notare abbiamo un solo valore per entrambe le variabili, dunque possiamo direttamente eliminare le due colonne e sostituire il nome della variabile “Value” con “Homicide_Rate_per100000”

omicidiWBT = omicidiWBT[,c(-3,-4)]
omicidiWBT <- rename(omicidiWBT, Homicide_Rate_per100000 = Value)
summary(omicidiWBT)

##  Country_Name       Country_Code            Year     
##  Length:15312       Length:15312       Min.   :1960  
##  Class :character   Class :character   1st Qu.:1974  
##  Mode  :character   Mode  :character   Median :1988  
##                                        Mean   :1988  
##                                        3rd Qu.:2003  
##                                        Max.   :2017  
##                                                      
##  Homicide_Rate_per100000
##  Min.   :  0.000        
##  1st Qu.:  1.487        
##  Median :  3.700        
##  Mean   :  8.418        
##  3rd Qu.:  9.500        
##  Max.   :139.132        
##  NA's   :12693

#Vincolo di chiave primaria
omicidiWBT %>% count(Country_Name, Year) %>% filter(n >1)

omicidiWBT %>% filter(is.na(Country_Name))

# Vincolo rispettato

Possiamo considerare ultimata la pulizia di questo dataset. Procediamo con l’importazione dei successivi.

2.2 Importazione file ottenuto da UNODC

Importiamo e puliamo il dataset su vari tipi di crimine ottenuto dal sito di UNODC, costituito da diversi fogli excel corrispondenti a diversi tipi di crimine

assault <- read_excel("CrimeUNODC.xlsx", sheet = 1, skip = 11)
glimpse(assault)

## Observations: 145
## Variables: 30
## $ Region              <chr> "Africa", NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ `Sub-region`        <chr> "Eastern Africa", NA, NA, NA, NA, NA, NA, ...
## $ `Country/territory` <chr> "Burundi", "Kenya", "Madagascar", "Mauriti...
## $ `2003`              <dbl> NA, NA, NA, 105, NA, NA, 19491, NA, 78464,...
## $ `2004`              <dbl> NA, NA, NA, 144, 1350, NA, 25744, NA, 9558...
## $ `2005`              <dbl> NA, 12715, NA, 122, 1272, NA, 35092, NA, 9...
## $ `2006`              <dbl> NA, 13186, NA, 106, 800, NA, 43782, NA, 77...
## $ `2007`              <dbl> NA, 12089, NA, 140, 724, NA, 28839, NA, 63...
## $ `2008`              <dbl> 437, 11479, NA, 148, 728, 2249, 21186, NA,...
## $ `2009`              <dbl> 496, 14189, NA, 274, 526, 1985, 17118, NA,...
## $ `2010`              <dbl> 304, 14078, 3450, 238, NA, 2004, 22979, NA...
## $ `2011`              <dbl> 499, 14366, 3434, 230, NA, 2526, 19023, NA...
## $ `2012`              <dbl> 455, 14534, 2115, NA, NA, 3177, 18030, NA,...
## $ `2013`              <dbl> 431, 13674, 2744, NA, NA, 3299, 16702, NA,...
## $ `2014`              <dbl> 549, 13949, 2542, NA, NA, NA, 14509, 1834,...
## $ `2015`              <dbl> NA, 14921, 2189, NA, NA, NA, NA, 2361, NA,...
## $ X__1                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ `2003__1`           <dbl> NA, NA, NA, 8.6835563, NA, NA, 74.3427356,...
## $ `2004__1`           <dbl> NA, NA, NA, 11.8437123, 6.5779347, NA, 94....
## $ `2005__1`           <dbl> NA, 35.969859, NA, 9.983584, 6.020824, NA,...
## $ `2006__1`           <dbl> NA, 36.33906892, NA, 8.63393266, 3.6802150...
## $ `2007__1`           <dbl> NA, 32.45322081, NA, 11.35444774, 3.237977...
## $ `2008__1`           <dbl> 4.9536404, 30.0148189, NA, 11.9546402, 3.1...
## $ `2009__1`           <dbl> 5.4280107, 36.1319183, NA, 22.0442591, 2.2...
## $ `2010__1`           <dbl> 3.2131513, 34.9084773, 16.3665873, 19.0712...
## $ `2011__1`           <dbl> 5.0969592, 34.6837662, 15.8403112, 18.3546...
## $ `2012__1`           <dbl> 4.494017, 34.163100, 9.486977, NA, NA, 29....
## $ `2013__1`           <dbl> 4.118113, 31.295716, 11.969697, NA, NA, 29...
## $ `2014__1`           <dbl> 5.075410, 31.092033, 10.784112, NA, NA, NA...
## $ `2015__1`           <dbl> NA, 32.401525, 9.032246, NA, NA, NA, NA, 4...

Notiamo che la struttura tidy non è rispettata. Inoltre abbiamo informazioni ridondanti in quanto le colonne da 2003 a 2015 rappresentano il numero di casi totali osservati nei diversi paesi, mentre le successive rappresentano i tassi. Con i prossimi comandi elimineremo le colonne superflue, rinomineremo quelle rimanenti e porteremo il dataset nella struttura tidy.

summary(assault$X__1)

##    Mode    NA's 
## logical     145

assault = assault[,-(4:17)]
glimpse(assault)

## Observations: 145
## Variables: 16
## $ Region              <chr> "Africa", NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ `Sub-region`        <chr> "Eastern Africa", NA, NA, NA, NA, NA, NA, ...
## $ `Country/territory` <chr> "Burundi", "Kenya", "Madagascar", "Mauriti...
## $ `2003__1`           <dbl> NA, NA, NA, 8.6835563, NA, NA, 74.3427356,...
## $ `2004__1`           <dbl> NA, NA, NA, 11.8437123, 6.5779347, NA, 94....
## $ `2005__1`           <dbl> NA, 35.969859, NA, 9.983584, 6.020824, NA,...
## $ `2006__1`           <dbl> NA, 36.33906892, NA, 8.63393266, 3.6802150...
## $ `2007__1`           <dbl> NA, 32.45322081, NA, 11.35444774, 3.237977...
## $ `2008__1`           <dbl> 4.9536404, 30.0148189, NA, 11.9546402, 3.1...
## $ `2009__1`           <dbl> 5.4280107, 36.1319183, NA, 22.0442591, 2.2...
## $ `2010__1`           <dbl> 3.2131513, 34.9084773, 16.3665873, 19.0712...
## $ `2011__1`           <dbl> 5.0969592, 34.6837662, 15.8403112, 18.3546...
## $ `2012__1`           <dbl> 4.494017, 34.163100, 9.486977, NA, NA, 29....
## $ `2013__1`           <dbl> 4.118113, 31.295716, 11.969697, NA, NA, 29...
## $ `2014__1`           <dbl> 5.075410, 31.092033, 10.784112, NA, NA, NA...
## $ `2015__1`           <dbl> NA, 32.401525, 9.032246, NA, NA, NA, NA, 4...

# riportiamo gli anni con la solita notazione eliminando i caratteri "__1" che compaiono alla fine.
colnames(assault) = gsub("__1","",colnames(assault))

# portiamo il dataset nella struttura tidy
assault = gather(assault, Year, Assault_rate_per100000, 4:16)

# convertiamo la colonna Year in formato intero
assault$Year = as.integer(assault$Year)
assault = rename(assault, Country_territory = `Country/territory`)
glimpse(assault)

## Observations: 1,885
## Variables: 5
## $ Region                 <chr> "Africa", NA, NA, NA, NA, NA, NA, NA, N...
## $ `Sub-region`           <chr> "Eastern Africa", NA, NA, NA, NA, NA, N...
## $ Country_territory      <chr> "Burundi", "Kenya", "Madagascar", "Maur...
## $ Year                   <int> 2003, 2003, 2003, 2003, 2003, 2003, 200...
## $ Assault_rate_per100000 <dbl> NA, NA, NA, 8.6835563, NA, NA, 74.34273...

# Verifica vincolo di chiave primaria
assault %>% count(Country_territory, Year) %>% filter(n > 1)

# Il vincolo di chiave primaria non è rispettato
assault = assault %>% filter(!is.na(Country_territory))

Abbiamo così importato i dati da un foglio excel che descrive il crimine “assault”. Riportiamo quindi di seguito le definizioni del suddetto crimine e degli altri tipi di crimine che verranno presi in considerazione:

1. Definitions: “Assault” means physical attack against the body of another person resulting in serious bodily injury, excluding indecent/sexual assault, threats and slapping/punching. ‘Assault’ leading to death should also be excluded.
2. Definitions: Total “Sexual violence” means rape and sexual assault, including Sexual Offences against Children.

Procediamo quindi con l’importazione dei rimanenti fogli excel di interesse.

sexViolence <- read_excel("CrimeUNODC.xlsx", sheet = 8, skip = 11)
glimpse(sexViolence)

## Observations: 132
## Variables: 30
## $ Region              <chr> "Africa", NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ `Sub-region`        <chr> "Eastern Africa", NA, NA, NA, NA, NA, NA, ...
## $ `Country/territory` <chr> "Burundi", "Kenya", "Madagascar", "Mauriti...
## $ `2003`              <dbl> NA, NA, NA, 298, NA, NA, NA, NA, NA, NA, N...
## $ `2004`              <dbl> NA, 2880, NA, 297, 834, NA, NA, NA, NA, NA...
## $ `2005`              <dbl> NA, 2590, NA, 314, 789, NA, 14101, NA, NA,...
## $ `2006`              <dbl> NA, 2965, NA, 420, 639, NA, 17046, NA, 0, ...
## $ `2007`              <dbl> NA, 3051, NA, 384, 629, NA, 18844, NA, 0, ...
## $ `2008`              <dbl> 599, 2719, NA, 413, 640, NA, 10365, NA, 0,...
## $ `2009`              <dbl> 700, 4068, NA, 442, 613, NA, 8646, NA, 0, ...
## $ `2010`              <dbl> 699, 4817, 696, 432, NA, NA, 8645, 9, 0, 2...
## $ `2011`              <dbl> 664, 4703, 711, 466, NA, NA, 8633, 0, 0, 2...
## $ `2012`              <dbl> NA, 4806, 596, NA, NA, NA, 9009, 6, NA, 19...
## $ `2013`              <dbl> NA, 4779, 419, NA, NA, 1734, 10974, 3, NA,...
## $ `2014`              <dbl> 1265, 5184, 893, NA, NA, NA, 13676, 22, NA...
## $ `2015`              <dbl> NA, 6164, 1115, NA, NA, NA, NA, 52, NA, NA...
## $ X__1                <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ `2003__1`           <dbl> NA, NA, NA, 24.644760, NA, NA, NA, NA, NA,...
## $ `2004__1`           <dbl> NA, 8.362986, NA, 24.427657, 4.063702, NA,...
## $ `2005__1`           <dbl> NA, 7.326932, NA, 25.695455, 3.734615, NA,...
## $ `2006__1`           <dbl> NA, 8.1711921, NA, 34.2099219, 2.9395718, ...
## $ `2007__1`           <dbl> NA, 8.1904853, NA, 31.1436281, 2.8131047, ...
## $ `2008__1`           <dbl> 6.7900014, 7.1095298, NA, 33.3599082, 2.78...
## $ `2009__1`           <dbl> 7.6604989, 10.3590559, NA, 35.5604471, 2.5...
## $ `2010__1`           <dbl> 7.38813398, 11.94446195, 3.30178108, 34.61...
## $ `2011__1`           <dbl> 6.7823264, 11.3544308, 3.2796917, 37.18810...
## $ `2012__1`           <dbl> NA, 11.29681143, 2.67339861, NA, NA, NA, 2...
## $ `2013__1`           <dbl> NA, 1.093771e+01, 1.827734e+00, NA, NA, 1....
## $ `2014__1`           <dbl> 11.6947062, 11.5550289, 3.7884391, NA, NA,...
## $ `2015__1`           <dbl> NA, 13.38536281, 4.60070995, NA, NA, NA, N...

sexViolence = sexViolence[,-(4:17)]
colnames(sexViolence) = gsub("__1","",colnames(sexViolence))
sexViolence = gather(sexViolence, Year, sexViolence_rate_per100000, 4:16)
sexViolence$Year = as.integer(sexViolence$Year)
sexViolence = rename(sexViolence, Country_territory = `Country/territory`)
summary(sexViolence)

##     Region           Sub-region        Country_territory       Year     
##  Length:1716        Length:1716        Length:1716        Min.   :2003  
##  Class :character   Class :character   Class :character   1st Qu.:2006  
##  Mode  :character   Mode  :character   Mode  :character   Median :2009  
##                                                           Mean   :2009  
##                                                           3rd Qu.:2012  
##                                                           Max.   :2015  
##                                                                         
##  sexViolence_rate_per100000
##  Min.   :  0.000           
##  1st Qu.:  6.287           
##  Median : 17.155           
##  Mean   : 33.856           
##  3rd Qu.: 46.624           
##  Max.   :419.788           
##  NA's   :631

# verifica vincolo chiave primaria
sexViolence %>% count(Country_territory, Year) %>% filter(n >1)

# eliminiamo i valori NA
sexViolence = sexViolence %>% filter(!is.na(Country_territory))

Passiamo ora alla pulizia e alla formattazione dei dati sull’immigrazione.

2.3 Importazione dati su immigrazione

immigrazione = read_excel("UN_MigrantStockTotal_2017.xlsx", sheet = 4, skip = 16, 
                          na = "..", col_names = FALSE)
head(immigrazione)

La situazione qui è più complicata: i nomi delle colonne nel foglio excel sono sfasati e di conseguenza si deve prcedere con un’assegnazione manuale dei nomi alle variabili guardando direttamente il foglio excel di interesse. Procediamo quindi con l’operazione:

colnames(immigrazione) = c("Sort_order","Area_Region_Country", "Notes", "Code",
                           "Type_of_data",
                           "1990", "1995", "2000", "2005", "2010", "2015", "2017",
                           "1990_m", "1995_m", "2000_m", "2005_m", "2010_m", "2015_m",
                           "2017_m", "1990_f", "1995_f", "2000_f", "2005_f", "2010_f",
                           "2015_f", "2017_f")
glimpse(immigrazione)

## Observations: 270
## Variables: 26
## $ Sort_order          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,...
## $ Area_Region_Country <chr> "WORLD", "More developed regions", "Less d...
## $ Notes               <chr> NA, "b", "c", "d", NA, "e", "e", "e", "e",...
## $ Code                <dbl> 900, 901, 902, 941, 934, 1503, 1517, 1502,...
## $ Type_of_data        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ `1990`              <dbl> 2.8614517, 7.1832325, 1.6766656, 2.1677033...
## $ `1995`              <dbl> 2.7940666, 7.8826379, 1.4927175, 1.9996841...
## $ `2000`              <dbl> 2.8088537, 8.6868948, 1.3964342, 1.5152279...
## $ `2005`              <dbl> 2.9123656, 9.6068724, 1.3923746, 1.2999352...
## $ `2010`              <dbl> 3.1620281, 10.5804404, 1.5609878, 1.177461...
## $ `2015`              <dbl> 3.3534532, 11.1913074, 1.7510442, 1.443059...
## $ `2017`              <dbl> 3.4133308, 11.5867310, 1.7762411, 1.440590...
## $ `1990_m`            <dbl> 2.8892949, 7.2359466, 1.7505268, 2.1845338...
## $ `1995_m`            <dbl> 2.8126583, 7.9299665, 1.5607685, 1.9988513...
## $ `2000_m`            <dbl> 2.8302973, 8.7381495, 1.4705551, 1.5198712...
## $ `2005_m`            <dbl> 2.9516978, 9.6640473, 1.4920500, 1.3268848...
## $ `2010_m`            <dbl> 3.2381085, 10.5476466, 1.7272761, 1.182836...
## $ `2015_m`            <dbl> 3.4267245, 11.0655808, 1.9293582, 1.442032...
## $ `2017_m`            <dbl> 3.4943529, 11.4607179, 1.9634866, 1.435144...
## $ `1990_f`            <dbl> 2.8341897, 7.1384502, 1.6008291, 2.1510868...
## $ `1995_f`            <dbl> 2.7761228, 7.8434658, 1.4228613, 2.0005757...
## $ `2000_f`            <dbl> 2.7880023, 8.6445955, 1.3203673, 1.5106542...
## $ `2005_f`            <dbl> 2.8733269, 9.5598606, 1.2898386, 1.2731602...
## $ `2010_f`            <dbl> 3.0856084, 10.6195078, 1.3896147, 1.172160...
## $ `2015_f`            <dbl> 3.2798463, 11.3190728, 1.5673068, 1.444112...
## $ `2017_f`            <dbl> 3.3318213, 11.7151751, 1.5833278, 1.446031...

Il dataset non è nella struttura tidy quindi è necessaria una riformattazione:

# creiamo la colonna Year e la variabile rappresentante la percentuale di migranti sulla popolazione totale
immigrazione = gather(immigrazione, Year, TotalPercentageOfMigrants, 6:12)
glimpse(immigrazione)

## Observations: 1,890
## Variables: 21
## $ Sort_order                <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1...
## $ Area_Region_Country       <chr> "WORLD", "More developed regions", "...
## $ Notes                     <chr> NA, "b", "c", "d", NA, "e", "e", "e"...
## $ Code                      <dbl> 900, 901, 902, 941, 934, 1503, 1517,...
## $ Type_of_data              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ `1990_m`                  <dbl> 2.8892949, 7.2359466, 1.7505268, 2.1...
## $ `1995_m`                  <dbl> 2.8126583, 7.9299665, 1.5607685, 1.9...
## $ `2000_m`                  <dbl> 2.8302973, 8.7381495, 1.4705551, 1.5...
## $ `2005_m`                  <dbl> 2.9516978, 9.6640473, 1.4920500, 1.3...
## $ `2010_m`                  <dbl> 3.2381085, 10.5476466, 1.7272761, 1....
## $ `2015_m`                  <dbl> 3.4267245, 11.0655808, 1.9293582, 1....
## $ `2017_m`                  <dbl> 3.4943529, 11.4607179, 1.9634866, 1....
## $ `1990_f`                  <dbl> 2.8341897, 7.1384502, 1.6008291, 2.1...
## $ `1995_f`                  <dbl> 2.7761228, 7.8434658, 1.4228613, 2.0...
## $ `2000_f`                  <dbl> 2.7880023, 8.6445955, 1.3203673, 1.5...
## $ `2005_f`                  <dbl> 2.8733269, 9.5598606, 1.2898386, 1.2...
## $ `2010_f`                  <dbl> 3.0856084, 10.6195078, 1.3896147, 1....
## $ `2015_f`                  <dbl> 3.2798463, 11.3190728, 1.5673068, 1....
## $ `2017_f`                  <dbl> 3.3318213, 11.7151751, 1.5833278, 1....
## $ Year                      <chr> "1990", "1990", "1990", "1990", "199...
## $ TotalPercentageOfMigrants <dbl> 2.8614517, 7.1832325, 1.6766656, 2.1...

# restano comunque le colonne relative alle percentuali maschili e femminili per ogni anno di osservazione. Risolviamo il problema ancora con la funzione gather:

immigrazione = gather(immigrazione, Year_m, MalePercentageOfMigrants, 6:12)
glimpse(immigrazione)

## Observations: 13,230
## Variables: 16
## $ Sort_order                <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1...
## $ Area_Region_Country       <chr> "WORLD", "More developed regions", "...
## $ Notes                     <chr> NA, "b", "c", "d", NA, "e", "e", "e"...
## $ Code                      <dbl> 900, 901, 902, 941, 934, 1503, 1517,...
## $ Type_of_data              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ `1990_f`                  <dbl> 2.8341897, 7.1384502, 1.6008291, 2.1...
## $ `1995_f`                  <dbl> 2.7761228, 7.8434658, 1.4228613, 2.0...
## $ `2000_f`                  <dbl> 2.7880023, 8.6445955, 1.3203673, 1.5...
## $ `2005_f`                  <dbl> 2.8733269, 9.5598606, 1.2898386, 1.2...
## $ `2010_f`                  <dbl> 3.0856084, 10.6195078, 1.3896147, 1....
## $ `2015_f`                  <dbl> 3.2798463, 11.3190728, 1.5673068, 1....
## $ `2017_f`                  <dbl> 3.3318213, 11.7151751, 1.5833278, 1....
## $ Year                      <chr> "1990", "1990", "1990", "1990", "199...
## $ TotalPercentageOfMigrants <dbl> 2.8614517, 7.1832325, 1.6766656, 2.1...
## $ Year_m                    <chr> "1990_m", "1990_m", "1990_m", "1990_...
## $ MalePercentageOfMigrants  <dbl> 2.8892949, 7.2359466, 1.7505268, 2.1...

immigrazione = immigrazione[,-15] #rimuovo colonna Year_m

immigrazione = gather(immigrazione, Year_f, FemalePercentageOfMigrants, 6:12)
glimpse(immigrazione)

## Observations: 92,610
## Variables: 10
## $ Sort_order                 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...
## $ Area_Region_Country        <chr> "WORLD", "More developed regions", ...
## $ Notes                      <chr> NA, "b", "c", "d", NA, "e", "e", "e...
## $ Code                       <dbl> 900, 901, 902, 941, 934, 1503, 1517...
## $ Type_of_data               <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Year                       <chr> "1990", "1990", "1990", "1990", "19...
## $ TotalPercentageOfMigrants  <dbl> 2.8614517, 7.1832325, 1.6766656, 2....
## $ MalePercentageOfMigrants   <dbl> 2.8892949, 7.2359466, 1.7505268, 2....
## $ Year_f                     <chr> "1990_f", "1990_f", "1990_f", "1990...
## $ FemalePercentageOfMigrants <dbl> 2.8341897, 7.1384502, 1.6008291, 2....

immigrazione = immigrazione[,-9] #rimuovo colonna Year_f

# vediamo la struttura
summary(immigrazione)

##    Sort_order    Area_Region_Country    Notes                Code       
##  Min.   :  1.0   Length:92610        Length:92610       Min.   :   4.0  
##  1st Qu.: 68.0   Class :character    Class :character   1st Qu.: 254.0  
##  Median :135.5   Mode  :character    Mode  :character   Median : 506.0  
##  Mean   :135.5                                          Mean   : 550.7  
##  3rd Qu.:203.0                                          3rd Qu.: 762.0  
##  Max.   :270.0                                          Max.   :5501.0  
##                                                                         
##  Type_of_data           Year           TotalPercentageOfMigrants
##  Length:92610       Length:92610       Min.   :  0.000          
##  Class :character   Class :character   1st Qu.:  1.476          
##  Mode  :character   Mode  :character   Median :  4.407          
##                                        Mean   : 11.560          
##                                        3rd Qu.: 13.905          
##                                        Max.   :100.000          
##                                        NA's   :735              
##  MalePercentageOfMigrants FemalePercentageOfMigrants
##  Min.   : 0.000           Min.   : 0.000            
##  1st Qu.: 1.351           1st Qu.: 1.293            
##  Median : 3.500           Median : 3.279            
##  Mean   : 8.953           Mean   : 8.368            
##  3rd Qu.:11.066           3rd Qu.:10.579            
##  Max.   :91.203           Max.   :87.858            
##  NA's   :11368            NA's   :11368

# convertiamo la colonna Year ad intero
immigrazione$Year = as.integer(immigrazione$Year)

Ora la struttura di questo dataset è adeguata: 9 variabili di cui le ultime 3 rappresentano le percentuali di “international migrant stock” sulla popolazione totale in aggregato e divisa per genere. Dove con “International migrant stock” si intende: “the total number of international migrants present in a given country at a particular point in time” (UN SD, 2017: 9, emphasis added). Andiamo a vedere se è rispettato il vincolo di chiave primaria. Come attributi della chiave consideriamo “Year” e “Code”

immigrazione %>% count(Year, Code)

# ci sono molti valori ripetuti
# consideriamo quindi i valori distinti
immigrazione = immigrazione %>% distinct(Year, Code, .keep_all = TRUE)

immigrazione %>% count(Year, Code) %>% filter(n > 1)

#ora il vincolo di chiave primaria è rispettato e non ci sono più tuple ripetute.

Per rendere informativa la colonna Type_of_data importiamo il foglio excel che contengono i valori delle note

note_immigrazione = read_excel("UN_MigrantStockTotal_2017.xlsx", sheet = 9, skip = 16, 
                         col_names = FALSE)

Fatto questo, prima di procedere con le analisi importiamo un dataset contenente le coordinate dei paesi, informzione che ci sarà utile per dei plot successivamente. La fonte del dataset è la seguente: https://developers.google.com/public-data/docs/canonical/countries_csv

country = read_delim("country.csv", delim = ";")
head(country)

# verifica vincolo di chiave primaria
country %>% count(country, name) %>% filter(n > 1)

# siccome il codice ISO rappresentato dalla variabile country sarà necessario per accedere alle coordinate rimuoviamo la tupla con codice ISO nullo
country = country %>% filter(!is.na(country))

3 ANALISI

Dopo aver importato, pulito e formttato i dati iniziamo con l’analisi. Cerchiamo quindi di ottenere un grafico che riporti l’andamento nel mondo, o almeno nei paesi principali, del fenomeno dell’immigrazione, per poi confrontare l’andamento della criminalità negli stessi stati.

3.1 Analisi geografica immigrazione

Iniziamo con un grafico semplice che mostri, per i diversi stati, l’andamento dell’immigrazione in funzione del tempo:

ggplot(immigrazione, aes(Year, TotalPercentageOfMigrants, group = Area_Region_Country)) +
  geom_line(alpha = 1/3) +
  labs(x = "Anno", y = "Percentuale di migranti") +
  theme_classic()

# molto poco indicativo

# cerchiamo di individuare solo i paesi soggetti ad immigrazione nel tempo
a = immigrazione %>%
  group_by(Area_Region_Country) %>%
  mutate(media = mean(TotalPercentageOfMigrants, na.rm = T)) %>%
  ungroup() %>%
  filter(media > 5 & media <20) %>%
  select(-media) %>% 
  ggplot(aes(Year, TotalPercentageOfMigrants)) + 
  geom_line(aes(col = Area_Region_Country), alpha = 0.4) +
  labs(x = "Anno", y = "Percentuale di migranti") +
  theme_classic()

ggplotly(a, width = 800, height = 600)

Per trovare i paesi maggiormente soggetti ad immigrazione si è deciso di considerare quelli che presentano una percentuale media negli anni compresa tra 5% e 20%. Questo perchè considerando valori medi troppo elevati risulterebbero molti paesi caratterizzati da una popolazione piccola, paesi per i quali è facile infatti che la percentuale media di migranti sulla popolazione totale sia elevata. Imponendo invece il limite inferiore di 5%, si escludono quei paesi che non sono stati soggetti ad un forte fenomeno di immgrazione negli ultimi tempi. Notiamo che non c’è nessun trend particolare per il fenomeno, molti paesi presentano un trend diverso e non si riconosce nessun particolare pattern. Con il grafico seguente invece andiamo a visualizzare la collocazione geografica e l’andamento dell’immigrazione per i paesi più soggetti a questo fenomeno, ma prima andiamo ad aggiungere una colonna che contenga il codice ISO a due caratteri, utile per i plot successivi, al dataframe “immigrazione”.

# aggiungiamo la colonna ISO
immigrazione = immigrazione %>% mutate(ISO = countrycode(Code, origin = "un",
                                                  destination = "iso2c"))
glimpse(immigrazione)

## Observations: 1,890
## Variables: 10
## $ Sort_order                 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...
## $ Area_Region_Country        <chr> "WORLD", "More developed regions", ...
## $ Notes                      <chr> NA, "b", "c", "d", NA, "e", "e", "e...
## $ Code                       <dbl> 900, 901, 902, 941, 934, 1503, 1517...
## $ Type_of_data               <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ Year                       <int> 1990, 1990, 1990, 1990, 1990, 1990,...
## $ TotalPercentageOfMigrants  <dbl> 2.8614517, 7.1832325, 1.6766656, 2....
## $ MalePercentageOfMigrants   <dbl> 2.8892949, 7.2359466, 1.7505268, 2....
## $ FemalePercentageOfMigrants <dbl> 2.8341897, 7.1384502, 1.6008291, 2....
## $ ISO                        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA,...

Procediamo quindi con i grafici successivi. Come già detto on il grafico seguente invece andiamo a visualizzare la collocazione geografica e l’andamento dell’immigrazione per i paesi più soggetti a questo fenomeno nel tempo .

a<- immigrazione %>%
  group_by(Area_Region_Country) %>%
  mutate(media = mean(TotalPercentageOfMigrants, na.rm = T)) %>%
  ungroup() %>%
  filter(media > 5 & media <20) %>%
  select(-media) %>%
  filter(!is.na(ISO)) %>% # L'ISO serve per accedere alle coordinate
  inner_join(country, by = c("ISO" = "country")) %>% 
  ggplot(aes(longitude, latitude, label=name)) +
  borders("world", alpha = 0.4) +
  geom_point(aes(size = TotalPercentageOfMigrants, frame = Year),shape = 21, alpha = 0.8, fill = "red") +
  coord_quickmap()+
  theme_classic()+
  labs(title="Immigrazione", x = "", y = "")
ggplotly(a, width = 900, height = 600) %>% animation_opts(frame = 1000, transition = 0, redraw = F)

Per questi paesi si nota un generale trend rialzistico e una collocazione non eccessivamente concentrata, il che dimostra la globalità del fenomeno. si può tuttavia vedere un’alta densità dei punti nella zona Europea. Ora invece andiamo a vedere come si è sviluppata l’immigrazione nel mondo, andando a vedere quanti paesi e quali hanno superato una certa soglia di immigrazione nel tempo (si mantiene comunque una soglia superiore per cercare di non considerare i paesi con popolazione estremamente bassa).

a = immigrazione %>% 
  filter(TotalPercentageOfMigrants > 7 & TotalPercentageOfMigrants <25) %>%
  filter(!is.na(ISO)) %>%
  left_join(country, by = c("ISO" = "country"))%>%
  ggplot(aes(longitude, latitude, label= name)) +
  borders("world", alpha = .4) +
  geom_point(aes(size = TotalPercentageOfMigrants, frame = Year), shape = 21, na.rm = T,      alpha =.8, fill = "red") +
  coord_quickmap()+
  theme_minimal()+
  labs(title="Immigrazione", x = "", y = "", size ="TotalPercentageOfMigrants")
ggplotly(a, width = 900, height = 600) %>% 
  animation_opts(frame = 500, transition = 0, redraw = F)

Si notano diverse zone particolarmente interessate dal fenomeno, e ad esempio in Europa si nota una tendenza col tempo verso un incremento di aree dei cerchi, indicndo dunque una maggior concentrazione di movimenti migratori.

3.2 Analisi geografica omicidi

Andiamo ora ad analizzare il dataset sugli omicidi preso da World Bank. Prima di tutto aggiungiamo una colonna contenente il codice identificativo ISO a due cifre tramite la funzione countrycode.

omicidiWBT =omicidiWBT %>% mutate(ISO = countrycode(Country_Code, origin = "iso3c", destination = "iso2c"))
glimpse(omicidiWBT)

## Observations: 15,312
## Variables: 5
## $ Country_Name            <chr> "Aruba", "Afghanistan", "Angola", "Alb...
## $ Country_Code            <chr> "ABW", "AFG", "AGO", "ALB", "AND", "AR...
## $ Year                    <int> 1960, 1960, 1960, 1960, 1960, 1960, 19...
## $ Homicide_Rate_per100000 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ ISO                     <chr> "AW", "AF", "AO", "AL", "AD", NA, "AE"...

Per questo dataset la serie storica inzia dal 1960 ma abbiamo molti valori NA per i primi decenni. Andiamo comunque a dare una visualizzazione grafica per gli stati con un tasso di omicidio più elevato nel tempo

# Usando dimensione dei punti
a = omicidiWBT %>%
  filter(Homicide_Rate_per100000 > 4) %>%
  filter(!is.na(ISO)) %>%
  inner_join(country, by = c("ISO" = "country")) %>%
  ggplot(aes(longitude, latitude, label= name, frame = Year)) +
  borders("world", alpha = .4) +
  geom_point(aes(size = Homicide_Rate_per100000), fill = "red",shape = 21, na.rm = T, alpha =.7) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="Tasso di omicidio")
ggplotly(a, width = 900, height = 600) %>% 
  animation_opts(frame = 500, transition = 0, redraw = F)

Vediamo che c’è un’alta concentrazione del fenomeno in diverse zone del globo, tra tutte spicca il centro-america, e a seguire abbiamo il sud-est asiatico e l’Africa (i cui dati compaiono da un certo istante in poi). Non mancano stati industrializzati e infatti vediamo che fra i punti compaiono gli Stati Uniti e alcuni stati della zona Europea.

Proviamo ora ad andare a graficare nella mappa mondiale i paesi con un alto tasso di omicidio e i paesi con un alto grado di immigrazione:

 b = immigrazione %>%
  filter(TotalPercentageOfMigrants > 7 & TotalPercentageOfMigrants <27, !is.na(ISO)) %>%
  inner_join(omicidiWBT %>%
               filter(Homicide_Rate_per100000 > 7, !is.na(ISO)), 
             by = c("ISO", "Year")) %>% inner_join(country, by = c("ISO" = "country")) 

a = ggplot(b, aes(frame = Year)) +
  geom_point(aes(longitude, latitude),col = "red", alpha = .6) +
  borders("world", alpha = .4) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="confronto omicidi e immigrazione", x = "", y = "", size ="")

ggplotly(a, width = 900, height = 600) %>%
  animation_opts(frame = 500, transition = 0, redraw = F)

Come possiamo notare non c’è una stretta corrispondenza tra immigrazione e tasso di omicidi in quanto solo pochi punti vengono rappresentati nel grafico, e quei punti rappresentano i paesi con alta percentuale di migranti e alto tasso di omicidio. Si nota una certa correlazione nella zona dell’america centrale e nelle isole della zona. Questa correlazione potrebbe essere dovuta a movimenti migratori verso stati più industrializzati per i quali i migranti sono costretti a fermarsi in stati caratterizzati da un elevata criminalità.

3.3 Analisi geografica crimine “aggressione”

Procediamo ora ad un’analisi dello stesso tipo per gli altri tipi di cirmine, cominciando dal crimine classificato come “Assault”, ovvero un attacco fisico verso un terzo che provoca ferite gravi, ad esclusione della morte e della violenza sessuale. Prima di tutto andiamo ad aggiungere al dataframe “assault” una colonna che contenga il codice ISO a due cifre per ottenere poi le coordinate dei paesi tramite un join con il dataframe “country”

assault = mutate(assault, ISO = countrycode(Country_territory, origin = "country.name", destination = "iso2c"))
# vediamo se ad ogni valore di Country_territory è stato assegnato uno ed un solo ISO
assault %>% count(ISO)

assault %>% count(ISO) %>% filter(n >13)

assault %>% filter(ISO == "GB")

Vediamo che il codice ISO “GB” si ripete per i diversi stati che compongono il Regno Unito, questo perchè viene assegnato lo stesso codice ISO a ciascuno stato che compone il Regno Unito. Tuttavia noi abbiamo bisogno che ci sia corrispondenza biunivoca tra nome dello stato e codice ISO, in modo da poter assegnare ad ogni stato una ed una sola coppia di coordinate geografiche. Per sistemare questo errore dobbiamo sostituire le 3 righe di ciascun anno con un unica riga che abbia come valore di “Assault_rate_per100000” la somma dei 3 valori delle righe da sostituire, essendo il tasso di aggressione del Regno Unito la somma dei tassi degli stati che lo formano.

# salviamo in una variabile di comdo un tibble che abbia come valore di  Assault_rate_per100000 di ciascun anno la somma dei 3 valori da rimuovere 
b1 =assault %>% group_by(Year, ISO) %>% filter(ISO == "GB") %>%  summarise(Assault_rate_per100000 = sum(Assault_rate_per100000))
# depuriamo il dataset originale dalle tuple con ISO pari a "GB" 
assault = assault %>% filter(ISO != "GB") 
# aggiungiamo alla variabile di comodo le colonne mancanti per poter unire le sue righe al dataset principale (assault)
b1 = b1 %>% mutate(Region = as.character(NA), `Sub-region` =as.character(NA), Country_territory = "United Kingdom")
# aggiungiamo le righe per ripristinare correttamente le tuple con ISO = "GB" precedentemente rimosse
assault = bind_rows(assault,b1)
assault %>% count(ISO) %>% filter(n >13)

Ora, con la solita metodologia, andiamo a vedere l’evoluzione di questo crimine nel mondo negli stati più soggetti a tale crimine. Andiamo quindi a graficare i paesi con un tasso di aggressione maggiore di una certa soglia nei rispettivi anni, in modo da identificare in quali paesi la situazione è più aggravata.

a = assault %>% filter(Assault_rate_per100000 > 65) %>% # 65 poco sopra la mediana
  inner_join(country, by = c("ISO" = "country")) %>%
  ggplot(aes(longitude, latitude, label= name)) +
  borders("world", alpha = .4) +
  geom_point(aes(size = Assault_rate_per100000, frame = Year), fill = "red",shape = 21, na.rm = T,   alpha =.8) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="Tasso di aggressione")

ggplotly(a, width = 900, height = 600)%>%animation_opts(transition = 0, redraw = F)

Anche in questo caso possiamo dire che il fenomeno è diffuso in molte regioni del globo con una maggior concentrazione in Amerca centrale, tuttavia si nota un calo negi ultimi anni. Bisogna però prestare attenzione perchè il calo osservato potrebbe esere dovuto ad un aumento di osservazioni NA. Detto questo, andiamo a vedere a livello grafico, come abbiamo fatto per il tasso di omicidi, se c’è una correlazione tra tasso di immigrazione e tasso di aggressione andando ad inserire nella mappa globale i punti corrispondenti agli stati che presentano contemporaneamente un alto tasso di aggressione ed un’alta percentuale di migranti.

b = immigrazione %>%
  filter(TotalPercentageOfMigrants > 7 & TotalPercentageOfMigrants <30, !is.na(ISO)) %>%
  inner_join(assault %>%
               filter(Assault_rate_per100000 > 65, !is.na(ISO)), 
             by = c("ISO", "Year")) %>% inner_join(country, by = c("ISO" = "country")) 

a = ggplot(b, aes(frame = Year, label = Area_Region_Country)) +
  geom_point(aes(longitude, latitude),col = "red", alpha = .6) +
  borders("world", alpha = .4) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="confronto aggressione e immigrazione", x = "", y = "", size ="")

ggplotly(a, width = 900, height = 600) %>%
  animation_opts(frame = 500, transition = 0, redraw = F)

Rispetto al crimine di omicidio, qui vediamo che c’è una correlazione grafica maggiore in quanto sono presenti diversi punti, in particolare in Europa.

3.4 Analisi geografica crimine “violenza sessuale”

Andiamo ora ad analizzare l’ultimo tipo di crimine: la violenza sessuale. Per il momento useremo la stessa metodologia utilizzata fino ad ora.

# inseriamo la colonna ISO
sexViolence = mutate(sexViolence, ISO = countrycode(Country_territory, origin = "country.name", destination = "iso2c"))

Per la stessa ragione di prima andremo a verificare se persiste il medesimo problema riscontrato per L’ISO del Regno Unito nei dataset assault e robbery.

# ci dovrebbe essere un solo codice ISO in corrispondenza ad ogni anno
sexViolence %>% count(Year, ISO) %>% filter(n>1)

Il problema è presente anche in questo dataset, andiamo dunque ad applicare la soluzione:

b1 =sexViolence %>% group_by(Year, ISO) %>% filter(ISO == "GB") %>% summarise(sexViolence_rate_per100000 = sum(sexViolence_rate_per100000))

sexViolence = sexViolence %>% filter(ISO != "GB") 

b1 = b1 %>% mutate(Region = as.character(NA), `Sub-region` =as.character(NA), Country_territory = "United Kingdom")

sexViolence = bind_rows(sexViolence,b1)

sexViolence %>% count(Year, ISO) %>% filter(n >1)

Come fatto in precedenza andiamo a visualizzare gli stati che presebtano una maggior concentrazione di questo crimine.

# violenze sessuali nel mondo mappando size
a = sexViolence %>% filter(sexViolence_rate_per100000 > 18) %>% # 18 è poco superiore alla mediana
  inner_join(country, by = c("ISO" = "country")) %>%
  ggplot(aes(longitude, latitude, label= name)) +
  borders("world", alpha = .4) +
  geom_point(aes(size = sexViolence_rate_per100000, frame = Year), fill = "red",shape = 21,
             na.rm = T, alpha = 0.8) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="Tasso di violenze sessuali", x = "", y = "", size ="Tasso di violenze sessuali")

ggplotly(a, width = 900, height = 600) %>% animation_opts(transition = 0, redraw = F)

A causa della mancanza di molti dati è difficile trarre qualche conclusione. Il numero di punti aumenta nel tempo ma in molti casi perchè per i primi anni non si era a conoscienza del valore dell’indice (NA nel dataframe). Guardando all’Europa la situazione non sembra chiara, all’inizio sembra che il fenomeno aumenti con una relativa diminuzione verso gli ultimi anni. Come per gli altri tipi di crimine andiamo a graficare nellamappa i paesi con più immigrazione e con più episodi di violenza sessuale.

# immigrazione eviolenze sessuali insieme
b = immigrazione %>%
  filter(TotalPercentageOfMigrants > 7 & TotalPercentageOfMigrants <27, !is.na(ISO)) %>%
  inner_join(sexViolence %>%
               filter(sexViolence_rate_per100000 > 18, !is.na(ISO)), 
             by = c("ISO", "Year")) %>% inner_join(country, by = c("ISO" = "country")) 

a = ggplot(b, aes(frame = Year, label = Area_Region_Country)) +
  geom_point(aes(longitude, latitude),col = "red", alpha = .6) +
  borders("world", alpha = .4) +
  coord_quickmap() +
  theme_minimal() +
  labs(title="confronto violenza sessuale e immigrazione", x = "", y = "", size ="")

ggplotly(a, width = 900, height = 600) %>%
  animation_opts(frame = 500, transition = 0, redraw = F)

Anche in questo caso solamente pochi stati risultano avere contemporaneamente un alto tasso di violenza sessuale e un’alta percentuale di immigrazione, si nota peraltro una concentrazione di punti non irrilevante in Europa.

4 REGRESSIONE

Quello che abbiamo fatto in precedenza è stata una semplice analisi grafica, con lo scopo di vedere quali stati erano maggiormente interessati dai fenomeni in esame. Con la parte seguente andremo a formalizzare l’assenza o la presena di relazione fra le quantità di interesse tramite un modello di regressione lineare, andando a verificare se effettivamente l’immigrazione spiega un aumento, o una diminuzione, della concentrazione dei crimini. Proviamo ora a fare quindi un’analisi più approfondita. Per fare questo eseguiamo un join fra il dataframe immigrazione e quelli riguardanti la criminalità consiederando solo gli stati in comune ad entrambi i dataframe (i valori NA verranno omessi).

4.1 Regressione tra immigrazione e tasso di omicidi

La regressione verrà sviluppata considerando come osservazioni la media del tasso di omicidi per ogni stato (media delle osservazioni annue), e la media della percentuale di migranti, in modo da avere per ogni stato una ed una sola osservazione. Procedendo in questo modo perdiamo molti punti, ma molti di questi sarebbero stati correlati fra loro (le osservazioni dello stesso stato sono simili) e avrebbero potuto portare ad una distorsione nell’elaborazione dei risultati.

# salviamo il dataframe su cui verranno effettuate le analisi
b = immigrazione %>% select(ISO, Year, TotalPercentageOfMigrants) %>%
  filter(!is.na(ISO)) %>%
  inner_join(omicidiWBT %>% filter(!is.na(ISO)), by = c("ISO", "Year")) %>%
  group_by(ISO) %>% summarise(mMigrants = mean(TotalPercentageOfMigrants, na.rim = T), 
                              mHomicide = mean(Homicide_Rate_per100000, na.rm = T)) %>%
  mutate(Country = countrycode(ISO, origin = "iso2c", destination = "country.name"),
         Continent = countrycode(ISO, origin = "iso2c", destination = "continent"))  

 a = b %>%  ggplot(aes(mMigrants, mHomicide)) +
  geom_point(aes(col = Continent, label = Country), alpha = .7, position = "jitter") +
  labs(x = "Tasso medio di migranti", y = "Tasso medio di omicidio") +
  theme_minimal()

ggplotly(a, width = 900, height = 600)

Questo grafico non risulta molto indicativo, c’è una massa di punti concentrata vicino all’origine e si possono notare molti valori estremi. Colorando i punti per continente si vede che si possono individuare dei gruppi all’interno dell’insieme di osservazioni. Proviamo con il seguente grafico a concentrare la visione sulla massa di punti vicino all’origine.

# limitando gli assi per una visione più chiara
a = a +
  xlim(0,50)+
  ylim(0,40) 

ggplotly(a, width = 900, height = 600)

Ora la situazione è più chiara, in ogni caso non è possibile notare una relazione netta fra le due grandezze. Vediamo che i paesi europei si concentrano molto vicino all’asse x, questo indica che in media negli ultimi anni sono stati caratterizzati da un’alta immigrazione ma da un basso tasso di omicidi. Vediamo invece che diversi paesi americani si avvicinano all’asse y, questo indica che tali paesi sono stati caratterizzati, in media nel corso degli ultimi anni, da relativamente poca immigrazione ma da un elevato tasso di omicidi. Proviamo ora a rimuovere i valori anomali o outliers.

# cercando di individuare i valori anomali per poi rimuoverli:
# creiamo una funzione che ci restituisca gli estremi di una distribuzione che identificano gli outliers
getOutliersBoundaries <- function(x){
  a = quantile(x, probs = 0.25, names = FALSE, na.rm = T) - 1.5*IQR(x, na.rm = T)
  b= quantile(x, probs = 0.75, names = FALSE, na.rm = T) + 1.5*IQR(x, na.rm = T)
  return(c(a,b))
}

#usiamo la funzione appena creata per filtrare i dati

b1 = b %>%
  filter(mMigrants > getOutliersBoundaries(mMigrants)[1] & mMigrants < getOutliersBoundaries(mMigrants)[2],
         mHomicide > getOutliersBoundaries(mHomicide)[1] & mHomicide < getOutliersBoundaries(mHomicide)[2]) 

a = b1 %>%  ggplot(aes(mMigrants, mHomicide)) +
  geom_point(aes(col = Continent, label = Country), alpha = .7, position = "jitter") +
  labs(x = "Tasso medio di migranti", y = "Tasso medio di omicidio") +
  theme_minimal()

ggplotly(a, width = 900, height = 600)

Dopo aver eliminato i valori anomali abbiamo un’immagine più chiara della situazione e vediamo una nuvola di punti che non sembra dare luogo ad una particolare relazione. Andiamo quindi a graficare un modello di regressione e a commentarne i parametri

# 1) Unica retta considerando anche gli outliers
a = b %>%
  ggplot(aes(mMigrants, mHomicide)) +
  geom_point(aes(col = Continent, label = Country), alpha = .5, position = "jitter") +
  labs(x = "Tasso medio di migranti", y = "Tasso medio di omicidio") +
  theme_minimal() 
  

a + geom_smooth(method = "lm")

# 2) una retta per continente considerando ancora gli outliers
a + geom_smooth(aes(col = Continent),method = "lm", se = F)

# 3) Unica retta rimuovendo outliers
a = b1 %>%
  ggplot(aes(mMigrants, mHomicide)) +
  geom_point(aes(col = Continent, label = Country), alpha = .7, position = "jitter") +
  labs(x = "Tasso medio di migranti", y = "Tasso medio di omicidio") +
  theme_minimal()

a + geom_smooth(method = "lm")

# 4) UNa retta per continente senza outliers
a + geom_smooth(aes(col = Continent),method = "lm", se = F)

Concentrandoci sui grafici che escludono gli ouliers vediamo che nel primo abbiamo una retta con pendenza negativa, ciò indicherebbe una correlazione tendenzialmente negativa fra tasso di omicidi e tasso di immigrazione. Guardando invece le rette differenziate per continente vediamo che l’Europa e l’Asia sono caratterizzate da rette praticamente orizzontali con una leggera tendenza positiva, le rette rappresentanti l’Oceania e l’America hanno una pendenza negativa, mentre l’unico contiente caratterizzato da una pendenza nettamente positiva è l’Africa. Andiamo ora ad analizzare alcuni parmetri del modello di regressione, prima considerando gli outliers e poi senza outliers

# Consiedrando gli outliers
mod1 = lm(mHomicide ~ mMigrants, data = b)
summary(mod1)

## 
## Call:
## lm(formula = mHomicide ~ mMigrants, data = b)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.704 -6.252 -3.233  1.302 74.227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.3638     0.9241  10.133   <2e-16 ***
## mMigrants    -0.1244     0.0510  -2.438   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.76 on 195 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.02959,    Adjusted R-squared:  0.02461 
## F-statistic: 5.945 on 1 and 195 DF,  p-value: 0.01565

# Non considerando gli otliers

mod2 = lm(mHomicide ~ mMigrants, data = b1)

summary(mod2)

## 
## Call:
## lm(formula = mHomicide ~ mMigrants, data = b1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.851 -3.872 -1.570  2.570 15.314 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.48809    0.51995  12.478   <2e-16 ***
## mMigrants   -0.10838    0.06406  -1.692   0.0926 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.731 on 159 degrees of freedom
## Multiple R-squared:  0.01768,    Adjusted R-squared:  0.01151 
## F-statistic: 2.862 on 1 and 159 DF,  p-value: 0.09263

Considerando gli outliers vediamo che il p-value per il test di nullità dei coefficienti è piuttosto basso, ciò indica che ci dovvrebbe essere relazione fra le due grandezze, e a giudicare dal segno del coeffiente relativo alla pendenza della retta di regressione, ci dovrebbe essere una relazione inversa. Peraltro vediamo che l’indice R^2 è molto basso, circa 3%, ciò indica che il modello non si adatta bene ai dati, infatti grazie a questo valore possiamo dire che appena il 3% della variabilità totale è spiegata dal modello. Eliminando invece gli outliers dalle osservzioni il coefficiente riguardante la pendenza della retta perde significatività e ha un valore leggermente negativo, indicando quindi una correlazione negativa fra immigrazione e tasso di omicidio. Anche qui tuttavia l’indice R^2 è basso, circa 1,7%, di conseguenza risulta difficile ammettere che valga una relazione lineare fra le due grandezze. Andiamo ora ad analizzare brevemente i residui dei due modelli.

# grafico residui modello con outlier compresi
b %>% 
   add_residuals(mod1) %>%
   ggplot(aes(mMigrants, resid)) +
   geom_point(aes(col = Continent), alpha = 0.3, position = "jitter") +
   geom_line(y = 0, col = "red") +
   theme_classic() +
   labs(title = "residui con outliers inclusi",x = "Tasso medio di migranti", y = "residui")

# Istogramma residui modello con outliers compresi
b %>%
   add_residuals(mod1) %>%
   ggplot(aes(resid)) +
   geom_histogram(binwidth = 3) +
   labs(title = "outliers inclusi", x = "Residui", y = "n") +
   theme_classic()

# grafico dei resiudi modello outliers esclusi
b1 %>% 
  add_residuals(mod2) %>%
  ggplot(aes(mMigrants, resid)) +
  geom_point(aes(col = Continent), alpha = 0.6) +
  geom_line(y = 0, col = "red") +
  theme_classic() +
  labs(title = "outliers esclusi", x = "Tasso medio di migranti", y = "residui", col = "Continente")

# Istogramma dei residui modello con outliers esclusi
b1 %>%
  add_residuals(mod2) %>%
  ggplot(aes(resid)) +
  geom_histogram(binwidth = 2) +
  labs(title = "outliers esclusi",x = "Residui", y = "n") +
  theme_classic()

Si nota che i residui sembrano confermare che non vi sia una relazione lneare fra le due grandezze. Dal primo istogramma si vede che la distribuzione non è normale ma presenta una lunga coda a destra, e ciò è confermato anche dallo scatter plot in cui si vede che la maggioranza dei punti giace al di sotto della retta di regressione stimata. Andando a rimuovere gli outliers la situazione migliora nello scatter plot, tuttavia dall’istogramma si può notare che la distrubuzione dei residui resta asimmetrica presentando un’assimetria positiva. Prima di proseguire proviamo ad utilizzare una regressione robusta per evitare di rimuovere gli ouliers senza però ottenere dei risultati distorti da questi ulitimi.

# definiamo il modello di regressione robusta
mod3 = MASS::rlm(mHomicide ~ mMigrants, data = b)
summary(mod3)

## 
## Call: rlm(formula = mHomicide ~ mMigrants, data = b)
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.292 -4.109 -1.559  3.501 76.667 
## 
## Coefficients:
##             Value   Std. Error t value
## (Intercept)  6.8985  0.4552    15.1561
## mMigrants   -0.0871  0.0251    -3.4652
## 
## Residual standard error: 5.806 on 195 degrees of freedom
##   (17 observations deleted due to missingness)

# confronto con il modello con gli outliers compresi
summary(mod1)

## 
## Call:
## lm(formula = mHomicide ~ mMigrants, data = b)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.704 -6.252 -3.233  1.302 74.227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.3638     0.9241  10.133   <2e-16 ***
## mMigrants    -0.1244     0.0510  -2.438   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.76 on 195 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.02959,    Adjusted R-squared:  0.02461 
## F-statistic: 5.945 on 1 and 195 DF,  p-value: 0.01565

# confronto grafico
b %>% 
  ggplot(aes(mMigrants, mHomicide)) +
  geom_point(alpha = 0.4) +
  geom_abline(aes(intercept = mod1$coefficients[[1]], slope = mod1$coefficients[[2]],
                  col = "ols")) +
  geom_abline(aes(intercept = mod3$coefficients[[1]], slope = mod3$coefficients[[2]],
                  col = "robust")) +
  labs(x = "Tasso medio di migranti", y = "tasso medio di omicidi") +
  scale_color_manual("Line colour", values = c("ols" = "red", "robust" = "blue")) +
  theme_classic()

Vediamo che tutto sommato i due modelli non si discostano molto nonostante ci sia una forte presenza di valori estremi. Si nota che nel modello robusto la pendenza della retta di regressione diminuisce leggermente accentuando il fatto che fra le due grandezze non ci sia una relazione lineare positiva. Per questo modello non si esegue la solita analisi dei residui e non si commenta l’indice \(R^{2}\) in quanto è costruito in maniera diversa dalla regressione semplice e non valgono le stesse proprietà. Siccome rimuovere gli outliers indiscrminatamente può rivelarsi una pratica pericolosa e non sempre corretta, siccome si sono notate rette di regressione diverse se differenziate per continente, e siccome nello scatter plot tasso medio di migranti vs tasso medio di omicidi si è notata una distribuzione dei punti con forma iperbolica proviamo nel seguito a trasformare i dati per verificare se effettivamente c’è una relazione iperbolica fra le due grandezze.

# grafico con ordinata trasformata (ordinata = tasso medio di omicidi)
a = b %>% 
  ggplot(aes(mMigrants, 1/mHomicide, col = Continent)) +
  geom_point(alpha = 0.5) +
  geom_smooth(aes(col = Continent), method = "lm", se = F) +
  labs(x = "tasso medio di migranti", y = "Reciproco tasso medio di immigrazione") +
  theme_light()
ggplotly(a, width = 900, height = 600)

# è necessario rimuovere i valori Nan e uguali a zero (con mHomicide > 0 si perdono 3 osservazioni) per sviluppare la regressione
c = b %>% filter(!is.nan(mHomicide), mHomicide > 0)

mod4 = lm(1/mHomicide ~ mMigrants + Continent, data = c)
summary(mod4)

## 
## Call:
## lm(formula = 1/mHomicide ~ mMigrants + Continent, data = c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09122 -0.18033 -0.04978  0.10021  2.22497 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.207756   0.055007   3.777 0.000213 ***
## mMigrants          0.011587   0.002052   5.647 5.94e-08 ***
## ContinentAmericas -0.211926   0.082818  -2.559 0.011286 *  
## ContinentAsia      0.229830   0.081545   2.818 0.005342 ** 
## ContinentEurope    0.331387   0.084330   3.930 0.000119 ***
## ContinentOceania   0.041684   0.126454   0.330 0.742039    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3931 on 188 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3436, Adjusted R-squared:  0.3261 
## F-statistic: 19.68 on 5 and 188 DF,  p-value: 9.516e-16

Vediamo che trasformando il tasso medio di omicidi nel suo reciproco la disposizione dei punti migliora, la relazione sembra puù lineare e questo indicherebbe una relazione iperbolica fra le grandezze originali. Si nota inoltre che tutti i coefficienti sono significativi tranne quello relativo al continente Oceania, ma soprattutto si vede che l’indice \(R^{2}\) è cresciuto notevolmente rispetto ai modelli precedenti. Andiamo allora ad analizzare brevemente i residui:

# scatter plot residui vs tasso medio di migranti
a = c %>% 
  add_residuals(mod4) %>%
  ggplot(aes(mMigrants, resid)) +
  geom_point(aes(col = Continent, label = Country), alpha = 0.5, position = "jitter") +
  geom_line(y = 0, col = "red") +
  theme_classic() +
  labs(x = "Tasso medio di migranti", y = "residui")
ggplotly(a, width = 900, height = 600)

# Istogramma dei residui
c %>%
  add_residuals(mod4) %>%
  ggplot(aes(resid)) +
  geom_histogram(aes(y =..density..), fill = "coral") +
  labs(x = "Residui", y = "density") +
  stat_function(fun=dnorm, color="red",
                args=list(mean=mean(mod4$residuals, na.rm = T), 
                          sd=sd(mod4$residuals, na.rm = T)))+
  theme_light() 

Si nota un vistoso miglioramento sia dallo scatter plot che dall’istogramma dei residui: la forma della distrubuzione è molto più simile ad una normale rispetto ai modelli precedenti, la variabilità è diminuita notevolmente, tuttavia rispetto ad una normle abbiamo una concentrazione molto più elevata in prossimità del valore zero. In conclusione possiamo dire che trasformando il tasso medio di omicidi nel suo reciproco e considerndo ne modello anche la variabile “Continente” otteniamo una relazione approssimativamente linere, in quanto la variabilità spiegata dal modello è circa pari al 34%, quindi un valore considerevole anche se tutto sommato non molto alto. Il fatto che la relazione sia abbastanza approssimativa si vede anche dal grafico dei resdui che di fatto non presenta le caratteristiche tipiche di un buon modello di regressione. Andiamo quindi a graficare il modello finale:

a = c %>% 
  add_predictions(mod4) %>%
  ggplot() +
  geom_point(aes(mMigrants, mHomicide, col = Continent), alpha = 0.5) +
  geom_line(aes(mMigrants, 1/pred, col = Continent)) +
  ylim(0,84) +
  labs(x = "tasso medio di migranti", y = "Reciproco tasso medio di immigrazione") +
  scale_colour_brewer(palette = "Set1") +
  theme_light()
ggplotly(a)

Fatto quetso andiamo ad analizzare la relazione fra l’immigrazione e gli altri tipi di crimine a partire con l’aggressione.

4.2 Regressione tra immigrazione e tasso di aggressione

Come prima andremo a fare un join tra i due dataframe di interesse per estrarre le informazioni che possono essere messe in comune.

b = immigrazione  %>%
  filter(!is.na(ISO)) %>%
  inner_join(assault, by = c("ISO", "Year")) %>%
  select(Area_Region_Country, Year, TotalPercentageOfMigrants, Assault_rate_per100000)

Come in precedenza analizzeremo l’eventuale relazione fra i tassi medi nel tempo per evitare di andare in contro a problemi di autocorrelzione fra i dati. Calcoliamo dunque i tassi medi:

b = b %>% group_by(Area_Region_Country) %>% 
  summarise(mMigrants = mean(TotalPercentageOfMigrants, na.rm = T),
            mAssault = mean(Assault_rate_per100000, na.rm = T)) %>%
   mutate(Continent = countrycode(Area_Region_Country, origin = "country.name", destination = "continent"))

Prima di tutto andiamo ad eseguire un plot di tutti i punti per vedere se visivamente c’è una relazione.

a = ggplot(b, aes(x = mMigrants, y = mAssault, col = Continent)) +
  geom_point(aes(label = Area_Region_Country), position = "jitter", alpha = .6, na.rm = T) +
  labs(x = "Percentuale media di migranti", y = "Tasso medio di aggressione") +
  geom_smooth(method = "lm", se = F, size = 0.5) +
  scale_colour_brewer(palette = "Set1") +
  theme_light() 

ggplotly(a, width = 900, height = 600)

Dal grafico si può vedere che le rette di regressione differenziate per continente sono fra loro diverse, questo ci suggerisce che la variabile Continent sarà significativa nel modello di regressione lienare. In ogni caso non sembra evidete nessuna chiara relazione lineare fra il tasso medio di migranti e il tasso medio di aggressione. Nel seguito andremo a commentare i modelli di regressione.

# modello di regressione lineare senza variabile Continent
mod1 = lm(mAssault ~ mMigrants, data = b)
summary(mod1)

## 
## Call:
## lm(formula = mAssault ~ mMigrants, data = b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -285.67 -166.54 -124.00   52.86 1540.58 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  174.390     31.420   5.550 1.59e-07 ***
## mMigrants      1.113      1.509   0.738    0.462    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 294.4 on 127 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.004266,   Adjusted R-squared:  -0.003574 
## F-statistic: 0.5441 on 1 and 127 DF,  p-value: 0.4621

# modello con variabile "Continent" compresa 
mod2 = lm(mAssault ~ mMigrants + Continent, data = b)
summary(mod2)

## 
## Call:
## lm(formula = mAssault ~ mMigrants + Continent, data = b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -345.70 -139.65  -90.36   27.75 1571.97 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        145.150     66.036   2.198   0.0298 *
## mMigrants            2.147      1.568   1.369   0.1734  
## ContinentAmericas  165.277     85.681   1.929   0.0560 .
## ContinentAsia      -57.748     84.778  -0.681   0.4970  
## ContinentEurope    -14.110     82.086  -0.172   0.8638  
## ContinentOceania    86.669    180.030   0.481   0.6311  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 287.6 on 123 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.08004,    Adjusted R-squared:  0.04265 
## F-statistic:  2.14 on 5 and 123 DF,  p-value: 0.06503

# modello con interazione
mod3 = lm(mAssault ~ mMigrants*Continent, data = b)
summary(mod3)

## 
## Call:
## lm(formula = mAssault ~ mMigrants * Continent, data = b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -438.32 -139.86  -96.23   76.56 1568.75 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                   -80.44     106.77  -0.753  0.45269   
## mMigrants                     120.98      45.09   2.683  0.00834 **
## ContinentAmericas             292.93     126.55   2.315  0.02234 * 
## ContinentAsia                 183.50     122.17   1.502  0.13576   
## ContinentEurope               223.85     119.86   1.868  0.06428 . 
## ContinentOceania              333.84     301.35   1.108  0.27017   
## mMigrants:ContinentAmericas  -101.77      45.71  -2.226  0.02787 * 
## mMigrants:ContinentAsia      -119.77      45.15  -2.653  0.00907 **
## mMigrants:ContinentEurope    -119.62      45.14  -2.650  0.00915 **
## mMigrants:ContinentOceania   -120.17      47.32  -2.539  0.01240 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 278.2 on 119 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.1671, Adjusted R-squared:  0.1041 
## F-statistic: 2.653 on 9 and 119 DF,  p-value: 0.007731

# test di significatività del modello
anova(mod2,mod3)

Dopo aver analizzato i parametri del modello di regressione possiamo dire che il modello migliore è quello che tiene conto della variabile Continent e dell’interazione che essa presenta con la percentuale media di migranti. Il modello ci dice che c’è una relazione significativa tra tasso medio di aggressione e tasso medio di migranti, tuttavia, a giudicare dal valore dell’indice \(R^{2}\) aggiustato pari appena a 10%, possiamo concludere che questo non si adatta bene ai dati, spiegando una bassa percentuale di variabilità. Con il grafico seguente andiamo a valutare visivamente l’inadeguatezza del modello andando a fare un grafico del tipo observed vs predicted

a = b %>% add_predictions(mod3) %>%
  ggplot(aes(x = mAssault, y = pred)) +
  geom_point(aes(col = Continent), alpha = 0.5) +
  labs(x = "osservati", "previsti") +
  geom_abline(intercept = 0, slope = 1, col = "red") +
  scale_colour_brewer(palette = "Set1") +
  theme_light()
  ggplotly(a, width = 900, height = 600)

Si vede che se il modello fosse stato corretto allora i punti si sarebbero disposti lungo la bisettrice del primo e terzo quadrante, qui la situazione non è per nulla simile a quella ideale. Questo conferma l’assenza di relazione lineare fra le quantità.

Proviamo anche questa volta a trasformare il tasso medio di aggressione nel suo inverso per vedere se si presentano miglioramenti.

# osservazioni che perderemo nella trasformazione
b %>% filter(mAssault == 0)

# scatter plot dopo la trasformazione
a = b %>% 
  ggplot(aes(mMigrants, 1/mAssault, col = Continent)) +
  geom_point(aes(label = Area_Region_Country), alpha = 0.5) +
  geom_smooth(aes(col = Continent), method = "lm", se = F, size = 0.5) +
  labs(x = "tasso medio di migranti", y = "Reciproco tasso medio di immigrazione") +
  scale_colour_brewer(palette = "Set1") +
  theme_light()
ggplotly(a)

Si nota che il grafico non risulta leggibile a causa di due valori anomali: l’Egitto e il Bangladesh. Questi due stati presentano un tasso medio di aggressione insolitamente basso vista la dimensione della popolazione dei paesi. A mio avviso sarebbe meglio rimuovere le osservzioni dal dataset visto il comportamento anomalo.

# elimino valori anomali
b = b %>% filter(Area_Region_Country != "Egypt",
                 Area_Region_Country != "Bangladesh")
# relazione lineare (?) per i diversi diversi continenti
a = b %>% 
  ggplot(aes(mMigrants, 1/mAssault, col = Continent, frame = Continent)) +
  geom_point(aes(label = Area_Region_Country), alpha = 0.5) +
  geom_smooth(aes(col = Continent), method = "lm", se = F, size = 0.5) +
  labs(x = "tasso medio di migranti", y = "Reciproco tasso medio di immigrazione") +
  scale_colour_brewer(palette = "Set1") +
  theme_light()
ggplotly(a, width = 900, height = 600) %>% animation_opts(frame = 1000, transition = 0, redraw = F)

Per alcuni continenti la relazione sembra essere approssimativamente lineare, andiamo dunque a costruire un modello di regressione.

# rimuovere valori nulli e Nan per costruire modello
c = b %>% filter(!is.nan(mAssault), mAssault > 0)
# regressione ordinata trasformata con variabile "Continent" inclusa
mod4 = lm(1/mAssault ~ mMigrants + Continent, data = c)
summary(mod4)

## 
## Call:
## lm(formula = 1/mAssault ~ mMigrants + Continent, data = c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11273 -0.03421 -0.00964  0.00706  0.62873 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.1104690  0.0252231   4.380 2.57e-05 ***
## mMigrants         -0.0007379  0.0006429  -1.148   0.2534    
## ContinentAmericas -0.0929809  0.0320250  -2.903   0.0044 ** 
## ContinentAsia      0.0051943  0.0322140   0.161   0.8722    
## ContinentEurope   -0.0679549  0.0307770  -2.208   0.0292 *  
## ContinentOceania  -0.0947993  0.0656718  -1.444   0.1515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1039 on 119 degrees of freedom
## Multiple R-squared:  0.1472, Adjusted R-squared:  0.1114 
## F-statistic: 4.108 on 5 and 119 DF,  p-value: 0.001778

# regressione ordinata trasformata con interazione
mod5 = lm(1/mAssault ~ mMigrants * Continent, data = c)
summary(mod5)

## 
## Call:
## lm(formula = 1/mAssault ~ mMigrants * Continent, data = c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12384 -0.03470 -0.00924  0.00719  0.62823 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                  0.13615    0.04355   3.126  0.00224 **
## mMigrants                   -0.01359    0.01764  -0.770  0.44271   
## ContinentAmericas           -0.11819    0.05059  -2.336  0.02121 * 
## ContinentAsia               -0.02113    0.04926  -0.429  0.66879   
## ContinentEurope             -0.09322    0.04892  -1.906  0.05918 . 
## ContinentOceania            -0.13221    0.11532  -1.147  0.25397   
## mMigrants:ContinentAmericas  0.01277    0.01787   0.715  0.47632   
## mMigrants:ContinentAsia      0.01289    0.01766   0.730  0.46704   
## mMigrants:ContinentEurope    0.01282    0.01768   0.725  0.46976   
## mMigrants:ContinentOceania   0.01358    0.01846   0.736  0.46353   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1054 on 115 degrees of freedom
## Multiple R-squared:  0.1513, Adjusted R-squared:  0.08484 
## F-statistic: 2.277 on 9 and 115 DF,  p-value: 0.02183

# test significatività di un modello
anova(mod4, mod5)

mod6 = lm(1/mAssault ~ mMigrants, data = c)
summary(mod6)

## 
## Call:
## lm(formula = 1/mAssault ~ mMigrants, data = c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.06150 -0.05294 -0.04052  0.00482  0.67343 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0630798  0.0122330   5.157 9.75e-07 ***
## mMigrants   -0.0004811  0.0006380  -0.754    0.452    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1104 on 123 degrees of freedom
## Multiple R-squared:  0.004602,   Adjusted R-squared:  -0.00349 
## F-statistic: 0.5687 on 1 and 123 DF,  p-value: 0.4522

I parametri ci dicono che non c’è relazione lineare fra il reciproco del tasso medio di aggressione e la percentuale media di migranti. Questo indica che la relazione tra tasso medio di aggressione e percentuale media di migranti non è nemmeno iperbolica. Dopo aver effettuato le seguenti analisi possiamo concludere che fra percentuale media di migranti e tasso medio di agressione non c’è relazione né iperbolica, né lineare.Inoltre la distribuzione dei punti non lascia intuire altre possibili relazioni. Andiamo quindi ad analizzare l’eventuale relazione fra migrazione e il crimine di violenza sessuale.

4.3 Regressione fra percentuale media di migranti e tasso di violenza sessuale

# join fra i dataframe di interesse
b = immigrazione  %>%
  filter(!is.na(ISO)) %>%
  inner_join(sexViolence, by = c("ISO", "Year")) %>%
  select(Area_Region_Country, Year, TotalPercentageOfMigrants, sexViolence_rate_per100000)

Useremo la stessa metodologia usata nei casi precedenti per evitare risultati distorti dall’autocorrelazione, ovvero andremo a computare la media delle grandezze d’interesse (in questo caso il numero di casi violenza sessuale su 100000 persone e la percentuale media di migranti) rispetto agli anni, in modo che ogni stato abbia una sola osservazione.

b = b %>% group_by(Area_Region_Country) %>% 
summarise(mMigrants = mean(TotalPercentageOfMigrants, na.rm = T),
mSexViol = mean(sexViolence_rate_per100000, na.rm = T)) %>%
mutate(Continent = countrycode(Area_Region_Country, origin = "country.name", destination = "continent"))

Andiamo ad eseguire il solito scatter plot delle grandezze di interesse con le rette di regressione differenziate per Continente. Se queste ultime presenternno un andamento simile allora si andrà a considerare una retta unica.

a = ggplot(b, aes(x = mMigrants, y = mSexViol, col = Continent)) +
  geom_point(aes(label = Area_Region_Country), 
             position = "jitter", alpha = .6, na.rm = T) +
  labs(x = "Percentuale media di migranti", y = "Tasso medio di violenza sessuale") +
  geom_smooth(method = "lm", se = F, size = 0.5) +
  scale_colour_brewer(palette = "Set1") +
  theme_light() 

ggplotly(a, width = 900, height = 600)

Le rette di regressione presentano una pendenza diversa, e in ogni caso anche in questo caso non c’è nessuna evidente relazione lineare fra le grandezze. Andiamo a vedere se questo è confermato dai parametri dei modelli di regressione.

# modello di regressione lineare senza variabile Continent
mod1 = lm(mSexViol ~ mMigrants, data = b)
summary(mod1)

## 
## Call:
## lm(formula = mSexViol ~ mMigrants, data = b)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -41.39 -27.83 -17.18  11.61 267.41 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.81857    4.94167   6.439 2.96e-09 ***
## mMigrants    0.09574    0.24356   0.393    0.695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44.31 on 114 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.001354,   Adjusted R-squared:  -0.007406 
## F-statistic: 0.1545 on 1 and 114 DF,  p-value: 0.695

# modello con variabile "Continent" compresa 
mod2 = lm(mSexViol ~ mMigrants + Continent, data = b)
summary(mod2)

## 
## Call:
## lm(formula = mSexViol ~ mMigrants + Continent, data = b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -58.004 -18.957  -9.059   3.210 268.751 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        10.4685    10.4507   1.002 0.318686    
## mMigrants           0.2191     0.2353   0.931 0.353716    
## ContinentAmericas  50.3099    13.0658   3.850 0.000198 ***
## ContinentAsia       0.8726    13.3234   0.065 0.947899    
## ContinentEurope    18.5791    12.5582   1.479 0.141881    
## ContinentOceania   50.5884    25.8185   1.959 0.052598 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 40.45 on 110 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.1968, Adjusted R-squared:  0.1603 
## F-statistic:  5.39 on 5 and 110 DF,  p-value: 0.0001798

# modello con interazione
mod3 = lm(mSexViol ~ mMigrants*Continent, data = b)
summary(mod3)

## 
## Call:
## lm(formula = mSexViol ~ mMigrants * Continent, data = b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.112 -18.335  -9.191   6.971 267.888 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                   -3.119     19.493  -0.160   0.8732  
## mMigrants                      9.804     11.769   0.833   0.4067  
## ContinentAmericas             44.893     21.724   2.067   0.0412 *
## ContinentAsia                 15.821     21.494   0.736   0.4633  
## ContinentEurope               36.313     20.953   1.733   0.0860 .
## ContinentOceania              44.034     44.103   0.998   0.3203  
## mMigrants:ContinentAmericas   -6.228     11.815  -0.527   0.5992  
## mMigrants:ContinentAsia       -9.669     11.773  -0.821   0.4133  
## mMigrants:ContinentEurope     -9.869     11.773  -0.838   0.4038  
## mMigrants:ContinentOceania    -8.335     11.940  -0.698   0.4866  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.05 on 106 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.2788, Adjusted R-squared:  0.2176 
## F-statistic: 4.553 on 9 and 106 DF,  p-value: 4.397e-05

# test di significatività del modello
anova(mod1, mod2)

anova(mod2,mod3)

Il modello migliore secondo i parametri di regressione è il terzo, ovvero quello che tiene conto dell’effetto della variabile Continent e dell’interazione. Vediamo che praticamente tutti i coefficienti non sono statisticamente significativi e questo andrbbe a confermare il fatto che fra le due grandezze non vi sia relazione lineare. L’indice \(R^{2}\) corretto è circa pari al 22%, dunque un valore basso ma tutto sommato non competamente trascurabile. Nel seguito andremo a fare un’analisi dei residui più approfondita rispetto ai modelli precedenti. Per poter iniziare costruiamo delle funzioni che ci saranno utili successivamente.

# indice dei 3 valori più grandi in un vettore
largestN = function(x, N) {
ndx <- order(x, decreasing = T)[1:N]  
return(ndx)
}
# indice dei 3 valori più piccoli in un vettore
smallestN = function(x, N) {
  ndx = order(x)[1:N]
  return(ndx)
}
# coefficiente angolare per costruire la qqline
qqslope = function(vec) {
  y <- quantile(vec[!is.na(vec)], c(0.25, 0.75))
  x <- qnorm(c(0.25, 0.75))
  slope <- diff(y)/diff(x)
  return(slope)
}
# intercetta della qqline
qqintercept = function(vec) {
  y <- quantile(vec[!is.na(vec)], c(0.25, 0.75))
  x <- qnorm(c(0.25, 0.75))
  int <- y[1L] - (diff(y)/diff(x)) * x[1L]
  return(int)
}

4.3.1 Diagnostica

# salviamo tutti i valori utili ottenuti dal modello di regressione lineare nella tibble d
d = as_tibble(augment(mod3)) %>% inner_join(b)
head(d)

# scatter plot residui vs tasso medio di migranti (linearità della relazione)
ggplot(d, aes(x = mMigrants, y = .resid, label = Area_Region_Country)) +
  geom_point(alpha = 0.4, col = "blue", position = "jitter") +
  geom_hline(yintercept = 0, linetype="dashed", col = "red") +
  geom_smooth(se = FALSE) +
  geom_text_repel(data = d[c(largestN(d$.resid, 3),smallestN(d$.resid,3)),],
                   nudge_y = 1) +
  labs(x= "percentuale media di migranti", y = "residui") +
  theme_light()

Vediamo che la disposizione dei residui non è ottimale, dovrebbero disporsi in maniera casuale e invece sembra che la maggior parte di questi giaccia sotto la retta \(x = 0\). Vediamo che i valori più estremi sono relativi agli stati: Regno Unito, Maldive, Svezia, Guatemala, Canada e Belize. Vediamo allora l’istogramma dei residui.

# Istogramma dei residui
d %>%
  ggplot(aes(.resid)) +
  geom_histogram(aes(y =..density..), fill = "coral") +
  labs(x = "Residui", y = "density") +
  stat_function(fun=dnorm, color="red",
                args=list(mean=mean(d$.resid, na.rm = T), 
                          sd=sd(d$.resid, na.rm = T)))+
  theme_light() 

Come già si intuiva dallo scatterplot la maggior parte dei residui è negativa, infatti si nota una forte concentrazione di punti al di sotto dello zero. Inoltre abbiamo l’occasione di notare una lunga coda a destra, sintomo di asimmetria positiva. Queste considerazioni mettono in dubbio la normalità dei residui. Di seguito allora utilizzeremo un altro strumento per valutare se i residui si distribuiscano effettivamente come una normale: il qqplot.

# qqplot (normalità della distribuzione)
a = ggplot(d, aes(label = Area_Region_Country)) +
  stat_qq(aes(sample = .std.resid), col = "blue", alpha = 0.5) +
  geom_abline(slope = qqslope(d$.std.resid),
              intercept = qqintercept(d$.std.resid), linetype="dashed", col = "red") +
  theme_light()

# procedimento per inserire delle etichette
d1<-ggplot_build(a)$data[[1]]
# si deve riordinare perchè ggplot_build ordina i valori
d1$Country = d$Area_Region_Country[order(d$.std.resid)]

ggplot(d1, aes(x =theoretical,y = sample, label = Country)) + 
  geom_point(alpha = 0.4, col = "blue") +
  geom_label_repel(data = d1[c(1,2,3,114,115,116),], nudge_y = 0.5) + 
  geom_abline(slope = qqslope(d$.std.resid),
              intercept = qqintercept(d$.std.resid), linetype="dashed", col = "red") + 
  theme_light()

Per costruire questo grafico dei residui si sono utilizzati i residui standardizzati, ovvero i residui trasformati in modo che abbiano media zero e vrianza unitaria. Anche da questo grafico è subito visibile il comportamento anomalo sulle code e vicino allo zero, tuttavia bisogna porre meno importanza alle code in quanto nel caso di un numero esiguo di dati, quindi come in questo caso, tendono ad essere instabili. Vediamo che abbiamo localizzato un ulteriore valore “strano”: Citta del Vaticano. Andiamo a vedere ora se nel nostro dataset sono presenti dei punti leva. Un punto si definisce leva in virtù del valore assunto dalle esplicative in quel punto rispetto agli altri, per effetto di queste, il modello tende ad avere un residuo (relativamente) piccolo in corrispondenza ad esso (ossia il modello è forzato a passarvi vicino). In sostanza un punto leva è un punto molto distante dalla massa di covariate, e già in precedenza abbiamo avuto modo di notare osservazioni di questo genere. Andiamo quindi ad utilizzare gli strumenti necessari per localizzarli.

# punti leva
# x vs h_ii
d %>% 
  ggplot(aes(mMigrants, .hat, label = Area_Region_Country)) +
  geom_segment(aes(xend = mMigrants, yend = 0), col = "red", size = 0.8, alpha = 0.4) +
  geom_label_repel(data = filter(d, .hat > 3*length(mod3$coefficients)/nrow(d)), nudge_y = 0.1) +
  geom_abline(intercept = 3*length(mod3$coefficients)/nrow(d), slope = 0) + #soglia critica
labs(x ="percentuale media di migranti", y = "h_ii") +
  theme_light()

Il grafico ci mostra diversi punti leva, alcuni erano già stati notati ma si sono aggiunti alla lista dei punti “particolari” anche “Bermuda” e “Emirati Arabi Uniti”. Andiamo a cercare di individuare i punti influenti, ovvero quei punti che determinano un grande cambiamento nelle stime dei parametri in caso di presenza/assenza dei punto stessi. Una misura di influenza è la distanza di Cook, tanto più graned è la distanza, tanto più influente sarà il punto.

# osservazioni influenti
# numero oss. vs distanza di Cook
d %>% 
ggplot(aes(as.numeric(.rownames), .cooksd, label = Area_Region_Country)) +
  geom_col(fill = "orange") +
  geom_label_repel(data = d[largestN(d$.cooksd,3),], nudge_y = 0.1) +
  labs(x ="numero osservazione", y = "Distanza di Cook") +
  theme_classic()

Abbiamo selezionato nel grafico i punti con distazna di Cook maggiore, vediamo che non sono nomi nuovi ma erano già stati nominati in precedenza per altre anomalie. Andiamo a commentare l’ultimo grafico per trovare i punti anomali, ovvero quelli che si discostano molto dal modello di regressione.

# osservazioni da evidenziare
d1= subset(d, .cooksd > .cooksd[order(-.cooksd)[4]] | .std.resid > 2.5 | 
             .hat > 3*length(mod3$coefficients)/nrow(d))
# leverage vs standard resid
ggplot(d, aes(.hat, .std.resid, label = Area_Region_Country)) +
  geom_point(aes(size=.cooksd), na.rm=TRUE, col = "blue", alpha = 0.5, show.legend = T) +
  stat_smooth(method="loess", na.rm=TRUE, col = "red", linetype="dashed", se = F,
              size = 0.5) +
  xlab("Leverage")+ylab("Standardized Residuals") +
  ggtitle("Residual vs Leverage Plot") +
  scale_size_continuous("Cook's Distance", range=c(1,5)) +
  geom_text_repel(data = d1, nudge_y = 0.1) +
  theme_light() +
  theme(legend.position="bottom")

Vediamo che in particolare in questo grafico andiamo a ritrovare tutti i punti che avevamo già commentato in precedenza. Facendo un riepilogo i punti più anomali sono: “United Kingdom” il quale si discosta molto dal modello e risulta essere anche influente, “Solomon Island” e “Holy See”, i quali risultano essere punti leva e influenti allo stesso tempo. Vediamo se eliminando queste osservzioni la qualità del modello cambia vistosamente.

b1 = b %>% filter(Area_Region_Country != "United Kingdom", 
                  Area_Region_Country !="Solomon Island",
                  Area_Region_Country != "Holy See")
mod3new = lm(mSexViol ~ mMigrants*Continent, data = b1)

tidy(mod3)

tidy(mod3new)

Vediamo che le considerazioni sul modello valgono anche dopo aver rimosso i valori critici, è confermata l’assenza di relazione lineare fra il tasso medio di violenze sessuali e la percentuale media di migranti.

# bontà del modello con valori anomali compresi
glance(mod3)

# bontà del modello senza valori anomali
glance(mod3new)

Si nota come rimuovendo i valori critici aumentà la bontà di adattamento del modello. Infatti si vede che l’indice \(R^2\) corretto è cresciuto vistosamente dopo aver rimosso le osservazioni particolari, raggiungendo un valore circa pari al 37%, quindi non affatto trascurabile.

5 Conclusione

In conclusione, una volta osservate le varie relazioni grafiche e i risultati ottenuti dalle regressioni, possiamo dire che non c’è una relazione diretta fra immigrazione e criminalità misurata in termini di tasso di aggressione e tasso di violenza sessuale. L’unica relazione che abbiamo colto è stata quella iperbolica fra tasso medio di omicidi e percentuale media di migranti che più che una relazione è un dato di fatto. Infatti sarebbe stupido dire che un basso numero di migranti causa un gran numero di omicidi, ma semplicemete gli stati più soggetti ad immigrazione sono quelli più industrializzati, peraltro caratterizzati da un basso tasso di omicidi. Al contrario, non sono mete dei movimenti migratori gli stati più poveri e instabili, putrtroppo caratterizzati da un elevato tasso di omicidi. Dunque la questione riguardante il rapporto tra immigrazione e criminalità rimane una questione delicata che richiede un grande impegno nella gestione del fenomeno, tuttavia possiamo concludere che questa questione è erroneamente oggetto di propoganda politica anti-immigrazione, in quanto non c’è significativa correlazione tra aumento degli stranieri e aumento dei crimini.